Modeling a circular equatorial test-particle in a Kerr spacetime. 
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Extreme Mass Ratio Inspirals (EMRIs) are one of the main gravitational wave (GW) sources for 
a future space detector, such as eLISA/NGO, and third generation ground-based detectors, like 
the Einstein Telescope. These systems present an interest both in astrophysics and fundamental 
physics. In order to make a high precision determination of their physical parameters, we need very 
accurate theoretical waveform models or templates. In the case of a circular equatorial orbit, the 
key stumbling block to the creation of these templates is the flux function of the GW. This function 
can be modeled either via very expensive numerical simulations, which then make the templates 
unusable for GW astronomy, or via some analytic approximation method such as a post-Newtonian 
approximation. This approximation is known to be asymptotically divergent and is only known up 
to 5.5PN order for the Schwarzschild case and to 4PN order for the Kerr case. A way to improve 
the convergence of the flux is to use re-summation methods. In this work we extend previous results 
using the Pade and Chebyshev approximations, first by taking into account the absorption of the 
GWs by the central black hole which was neglected in previous studies, and secondly by using the 
information from the Schwarzschild and absorption terms to create a Kerr flux up to 5.5PN order. 
We found that these two additions both improve the convergence. We also demonstrate that the 
best re-summation method for improving the flux model is based on a flux function which we call 
the "inverted Chebyshev approximation". 

I. INTRODUCTION. 

The Extreme Mass Ratio Inspirals (EMRIs) are one of the most important sources of gravitational waves (GWs) for 
future space detectors such as the eLISA/NGO mission [HQ, and terrestrial detectors like the Einstein Telescope 
These systems are composed of a stellar mass compact object orbiting around a (super)massive black hole (BH) such 
that the mass ratio is on the order of 10 -5 — 1CP 7 . Their physical interest comes from the fact that the small object 
normally traces out approximately 10 5 orbits from when it is first observed in the detector until it plunges into the 
central BH. The test particle can then be thought of as something which maps the spacetime around the central BH, 
providing a very strong test of the theory of General Relativity. 

Even if EMRIs are quite simple objects, their dynamics, which we need to understand for GW astrophysics, are 
very complicated. As we do not have a solution for the case of generic orbits around the central BH, we concentrate 
on a test particle in a circular equatorial orbit. Due to the fact that there is no plane precession in this case, we have 
one less constant of motion and the trajectories are completely described by the energy and angular momentum of 
the system. One of the key aspects in understanding the dynamics is the flux of GWs emitted by the system. At 
present, we have an exact expression for the binding energy of the system, but for the flux we use a post-Newtonian 
(PN) approximation. This approximation is known it up to 4PN for the Kerr case 043 and up to 5.5PN in the 
Schwarzschild case [§-fi~2| . 

While it was shown that the PN series breaks down for orbital separations of r < 10 M [HI, the vali dity of the 
PN flux was studied in the non-spinning and spinning cases respectively using an asymptotic analysis Il5j . These 
studies found that for the Schwarzschild case the best PN expansion is at 3PN, but a conclusion could not be made 
for the Kerr case due to the lack of higher order terms in the series. The second result from these studies was that 
the PN approximation is more suitable for retrograde than for prograde orbits. This is most easily explained by the 
position of the last stable orbit (LSO). When the dimensionless spin parameter, q, is positive, the LSO moves closer 
to the the central black hole as q increases, reaching a maximum value of r = M for 5=1, and so, the system is more 
relativistic causing the PN approximation to break down sooner. The contrary happens for q < as the position 
of the LSO moves outwards from the BH to a maximum value of r = 9M at a value of q = — 1. And because the 
PN approximation is an expansion in the invariant PN velocity parameter, x, and assumes that x << c, it makes 
sense that it is more convergent for less relativistic systems. Another interesting result is that the PN flux is so bad 
at small separations between the two bodies that it should not be used for r/M smaller than a limit in the interval 
r/M £ [6,54] depending on the spin of the central BH. These results, along with previous efforts to improve the 
convergence of the PN series, are our motivation to resolve this problem. 

Given the starting information provided by the PN approximation, a way to improve the convergence of a series 
is to use re-summation methods. The first proposed solution for the GW flux was to use a Pade approximation, 
which replaces the PN series by a rational polynomial. This was applied first to Schwarzschild orbits [16|, and was 
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later extended to Kerr orbits [17]. In both cases, they found that in general, the Pade flux has a better and more 
monotonous convergence than the PN flux. Even so, the Pade solution was not complete due to the existence of poles 
in the approximation at certain PN orders, and also due to the fact that at some PN orders the new approximation 
was worse than the original PN series. In previous articles it was demonstrated that one could create both sub and 
super-diagnonal Pade approximations, depending on the number of terms in either the numerator or the denominator. 
The problem of poles could be cured by moving to a sub-diagonal approximation should a pole occur in the super- 
diagonal approximation and vice versa. However, due to the inherent form of the approximation, there is no cure 
for the cases where the approximation simply performs badly. Furthermore, in the Kerr case, a variety of spins were 
investigated, but the flux was kept at the 4PN order, which we now know to be divergent. The second solution which 
was proposed was to replace the divergent power series by a Chebyshev series, and was tested in the non-spinning 
test-mass and comparable mass cases 

[laliai- As a Chebyshev series is closely related to the minimax polynomial, 
they were demonstrated to obtain a better convergence than using the PN or Pade fluxes. This improvement was due 
to two aspects : the first is that while power series have a radius of convergence that may not cover the entire domain 
of interest, Chebyshev series have an ellipse of convergence that spans the entire domain of interest. Secondly, as we 
shall see later, the PN factorized flux contains logarithmic terms which begin at the 3PN order. This part of the 
series is untreatable with Pade approximation as it yields coefficients that are either null or infinite. The Chebyshev 
approximation has no problems with the fact the the first non-zero term of the series is at 3PN order and we can thus 
define a Chebyshev series here as well. 

In this work we extend these previous results, first, by taking into account the absorption of the GW by the central 
black hole. Secondly, in the spinning case, as the 4PN flux is divergent, we include the information that we have in 
the Schwarzschild case up 5.5PN in the Kerr flux function. Then, we search for the best factorized flux function, from 
the function proposed in Ref. [l6| where the velocity at the LSO is used as a normalising velocity, and of some of 
its variants, which include using the velocity at the photon ring as a normalising velocity or inverting the factorized 
flux function. Finally, we study the convergence of this more complete PN flux, and of the Pade and Chebyshev flux 
models, for different orders of approximation and different spins. 

We will see that the flux needs to be treated differently for different types of orbits : the non-spinning, retrograde and 
prograde cases. This separation stems from the differing relativistic nature involved as the test-particle approaches 
the central black hole and will influence our choice of normalizing velocity in the flux function. We will use the 
invariant velocity at the photon ring for the retrograde case, and the invariant velocity at the last stable orbit for the 
non-spinning and prograde cases. We will also see that in most situations, it is better to use an inverted factorized 
flux, introduced in Ref 17 1 , rather than the original factorized flux, proposed by Ref [lj| . And finally we will see that 
for the Chebyshev re-summation, we will not necessarily use the highest known PN order, especially if that order less 
convergent than at lower orders. 

This work is organized as follow, in Section [TT| we will describe the PN flux. In Section IIII1 we will describe the 
Pade and Chebyshev approximations, and we will show how to construct these respective fluxes. In Section IIV1 we 
will present our results for the different spins. Finally, in Section [V] we will provide our conclusion. 



II. POST-NEWTONIAN GW FLUX. 

Search algorithms for GW astronomy are based on the technique of matched filtering [2(| , which cross-correlates 
the received signal with theoretical templates of known parameter values. In order to determinate as precisely as 
possible the physical parameters of the observed source, we need very accurate templates. As matched filtering is 
very sensitive to the GW phase, it is therefore imperative that we have a very accurate phase model. The GW phase 
can be written in the form 

</>(«) =0o + 2/ v 3 —j—t-dv, (1) 
Jv F(v) 

where <j)Q is a constant reference GW phase (which can be taken to be when the signal is first seen in the detector), 
E(v) is the binding energy and E'(v) = dE(v)/dv. F(v) is the flux of the GW, and v is the local linear velocity given 
by v = \J M/r, where r is the radial coordinate in Boyer-Lindquist coordinates and M is the mass of the central black 
hole. In the case of equatorial and circular orbits, we have an exact analytical expression for the binding energy pH | 

VI - 3v 2 + 2qv 3 

where r\ = /j,M/(/j, + M) 2 ~ [ijM is the symmetric mass ratio, /i is the mass of the stellar mass object, and q is the 
dimensionless spin parameter defined as q = a/M. 
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From this expression we can determine three particular orbits. First, we can determinate the position of the last 
stable orbit (LSO) which corresponds to the condition E'(y) = [2l[ 



rL(l)= M 3 + z 2 (q) T x/[3 - zM] [3 + zi(<z) + 2 z 2 (q)} 



where 



Zl(q) = l + (l-q 2 V (l + q^ + (l-q) 



(3) 



(4) 



22(g) = \/3g 2 



(5) 



The signs + and — correspond to prograde and the retrograde orbits respectively. Then we can obtain the last 
unstable orbit, also known as the photon ring, where the derivative of the binding energy, E'(v), has a pole [2lL l22l| ■ 



r%{q)=2M 



1 + cos 



cos 1 (=Fg) 



Finally, we can determine the position of the black hole horizon, where the energy is zero 



r% =M (1 + VT 



(6) 



(7) 



The position of these three orbits varies with the sign and magnitude of q. In the Schwarzschild case the horizon, 
the photon ring and the last stable orbit are respectively at 2A/, 3M and 6M. When q is positive and increasing, 
these orbits move to smaller radius until the coincide at radius r = M in the limiting case of a maximally spinning 
BH, i.e. q = 1. When q is negative and decreasing, the LSO moves away from the black hole, while the photon ring 
and the horizon again approach the central BH. 

While we have a exact expression for the binding energy in the case of circular equatorial orbit, it is not the case 
for the GW flux function F(v). Instead we use a PN approximation where we have a series approximation in powers 
of v/c. This approximation has been calculated up to 5.5PN order in the case of a Schwarzschild BH and up to 4PN 
in the case of a Kerr BH. The general form of the flux in both cases is given by the expression [5I-H2I] 



F Tri {x-q) =F N (x) 



^ ak{q)x k + \n(x) ^2 M?) 5 



fc=0 



fe=6 



where n — 8 and 11 in the Kerr and Schwarzschild cases respectively, and Fjy(x) is the Newtonian flux given by 

F N (x) = B rfx™. 
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(8) 



(9) 



Here, x is the invariant PN velocity parameter observed at infinity, which is related to the orbital angular frequency 
f2 by x = (MfJ) 1 / 3 (Note that this is a different use of notation to the comparable mass literature where we define 



(Mil) 2 ^ 3 ). The relation between the invariant velocity x and local linear velocity v is given by [23 

x(v, q) = v(l - qv 3 + q 2 v 6 ) 1/3 , 



(10) 



which reduces to x — v in the Schwarzschild limit. The coefficients au{q) and bk(q) appearing in the flux function 
are given in Appendix [X] They are the sum of three contributions : a spinless term, a spin dependent term, and an 
absorption term (23|. The absorption term comes form the absorption of the GW by the central BH. 

Until now, in the Kerr case, the two series in Eqn ([5]) have been stopped at n = 8 as we do not know the spin- 
dependent term for higher orders. However, this 4PN flux function is divergent. In this work we propose an extension 
of this series up to n = 11 where we use the spin-independent and the absorption terms for n > 8. While this flux 
function is still incomplete, due to the lack of knowledge of spin dependent coefficients at greater than 4PN order, 
we shall demonstrate that this flux function has a better convergence when compared with a numerical flux that the 
truncated 4PN Kerr flux function. 
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III. APPROXIMATING THE GW FLUX FUNCTION. 



The field of numerical analysis presents a number of ways of improving the convergence of a series. Pade approx- 
imation is a powerful method due to the fact that it is the best way of representing a power series by means of a 
rational polynomial, i.e. 



N 



T n {x) =J2 a » xk » *W = (11) 

d kx k 



k=0 

fe=0 



with n + l = N + M + l and the power series of the right hand side being identical to the original power series. 
Theory suggests that if we think of the Pade approximation as a matrix, the best convergence is achieved by staying 
as close to the diagonal as possible, i.e. keep the number of terms in the numerator and denominator equal. At odd 
orders of approximation, it is possible to allow either the number of terms in the numerator (super- diagonal) or the 
denominator (sub-diagonal) increase by one term. Thus, if we choose M = N + e, where e = or 1, we have the 
diagonal and sub-diagonal approximations, respectively. The advantage of the sub-diagnoal Pade approximant is that 
we can write it in the form of a continued fraction 

N 

pN _ i=0 _ CO / 19 ^ 



C2X 



1 + C n X 



In this form, when we go from the order n to the order n+1, only the (n+1)*' 1 coefficient needs to be calculated [24.l25j. 
all lower order coefficients remain constant. In contrast, with a super-diagonal approximation, each higher order of 
approximation requires a re-calculation of all coefficients. Pade approximation is a very useful method for accelerating 
the convergence of a series as not only can it well represent zeros and poles, but it can provide new information 
regarding the underlying function. The main problem with Pade approximation is when it gets the approximation 
wrong. As with any power series, the Pade approximation is ultimately based on the Weierstrass theorem which only 
guarantees convergence with a high number of terms, it is impossible to know at any particular order if the Pade 
approximation will perform well. 

On the other hand, the best way of increasing the convergence of a power series by using a polynomial representation 
is to use the minimax polynomial, which unfortunately is not that easy to find. The minimax polynomial approximates 
the underlying function such that we are presented with an oscillating equal error curve. This approximation has the 
ability to allow the error to grow in parts of the domain where we already have high accuracy in order to improve 
the accuracy in parts of the domain where the error begins to grow. Closely related to the minimax polynomial is a 
polynomial based on Chebyshev polynomials of the first kind which can be defined using the iterative relation 



with 



T n (x) = 2xT n ^ 1 (x) - T n - 2 (x), (13) 



T (x) = l , T l {x)=x. (14) 



The Chebyshev polynomials are orthogonal polynomials and are a special case of the Gegenbauer polynomials (see 
Ref [HI, Hg] for a more in depth explanation). They are defined in the interval [—1,1], where at the n th order they 
have n zeros and n + 1 extrema such that |T n (x)| = 1. As most problems are not limited to this interval and are 
defined in the interval [a, b], so we define the shifted Chebyshev polynomials of the first kind, T*(x) = T(s), by 

s= 2x-(a + b) Vie[o,i],s6[-l,l]. (15) 
o — a 

The shifted Chebyshev polynomials can be calculated using the following iterative expression 

T:(x) = 2 ( 2X ~ {a + b) ) T^x) - T:_ 2 (x), (16) 
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with 

!?(*)= 1 , 3?(*) = ^^. (17) 

The Chebyshev approximation again tries to minimize the maximum error, albeit less successfully than the minimax 
polynomial, but sufficiently so to give a more convergent series at lower orders of approximation than the truncated 
power series. While a Pade approximation can provide new information on an underlying function, a Chebyshev 
approximation essentially re-packages the information that we have from the power series into a more useful form. 
Another beauty of the Chebyshev series is that each coefficient at higher order is smaller than the previous coefficient. 
As the Chebyshev polynomials have a maximum magnitude of unity, this means that the truncation error associated 
with a Chebyshev series is essentially the value of the truncated coefficient. Ultimately, the performance of the 
Chebyshev series is dependent on the generating power series. If a series at a particular order is particularly divergent, 
then the corresponding Chebyshev series with perform better, but will still be divergent. We will see that in this 
case, starting at a lower order of approximation in the power series can result in a more convergent Chebyshev 
approximation. 



A. Pade GW flux approximation. 

The first method that was introduced in order to try and improve the convergence of the flux function was based 
on a Pade approximation [l6| . Due to a number of issues with the PN approximation to the flux, a direct application 
of Pade approximation is not possible. The first problem with the PN flux are the logarithmic terms which appears 
at 3PN, which effect the convergence of the approximation. In order to reduce this effect, it was shown that one 
could normalize the velocity parameter appearing in the logarithm and also factor out the logarithmic terms from the 
original series. This led to the creation of a new flux function 



F Tn (x) = F N (x) 



where 



1 + ln 



x N 



fe=6 



J2c k x k , 



(18) 



. fc=0 



Cfc 



a-k 



\n(x N )b k 



V k < 6 

V k > 6 



(19) 



and 



fc-6 



(20) 



In Ref [l6j], the invariant velocity at the last stable orbit is used as the normalizing velocity (i.e. xn — %lso)- This 
value was used as the intention was to model the inspiral of the test particle only as far as the LSO. As we approach 
this point the effect of the logarithmic terms are greatly diminished. However, as there are no constraints in choosing 
the value of the normalizing velocity, we will also investigate other possibilities. 

A second issue is that the linear term of the PN flux is null, which creates a problem when we use the Pade 
approximation. A null term in a series can result in null and infinite Pade coefficients when we write the expansion 
in continued fraction form. Furthermore, we know that the flux function has a pole at the photon ring, so in order 
that we use the full potential of the Pade approximation it was required that the function exhibit this pole [24|, H3] ■ 

Therefore, by multiplying the flux by ^1 — jr-^ we obtain the factorized flux 



fr N = F N (x) 



where 



— >D 

fe=6 



kX 



(21) 



J fc=0 



/o 

fk 



Co, 
Cfc - 



' pr 



(22) 
(23) 
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and where x pr is the invariant velocity at the photon ring. 

There exits another factorized flux introduced in Ref [171 ] ■ that we will call the inverted factorized flux, which is 
obtained by analytically inverting the previous equation, which gives, 



flTjx) 



F N (x) 



1 I — li-r L 

XUo ' k=6 



k=0 



(24) 



where 



fc=0 



\k=0 



and the dfc coefficients are calculated using 



d 
dk 



1 

To' 



ft 



1 



(25) 

(26) 
(27) 



As the sum over the Z^'s in Eqn (|2"Tj) and Eqn (IM1) only starts at fc = 6, we cannot apply the Pade approximation as 
it would generate null, undefined and infinite coefficients in the Pade approximation. So the Pade approximation is 
only applied to the non-logarithmic sum in the above equations, and we hence obtain the Pade flux [16| 



FpJx) =F N (x) ( 1 



1 + h JL 

v 7 k=6 



pN 



E/< 

Lfc=0 



(28) 



and the inverted Pade flux \l'< 



F IPn (x) = F N (x) 



1 - 



i-M— )E** a 



PS 



n-1 



,k=0 



d k x k 



(29) 



The Pade approximant was studied in the non-spinning [l6[ and spinning [l?} cases, and it was shown that in both 
case the Pade flux better modeled a numerical flux than the PN flux. However it was also seen in both cases that 
the Pade approximation can suffer from singularities at different PN orders. In principle we can recover from this by 
moving to a super-diagonal approximant. Furthermore, in some cases the Pade approximation was a worse fit than 
the original PN series. This was due to the fact that the factorized flux is highly divergent at certain PN orders and 
performed much worse than the PN series. While the Pade approximation was able to improve the convergence of 
the factorized flux, it still remained less convergent than the PN flux. This also turned out to be the case in the Kerr 
study. The convergence of the flux was studied for different spin values of the central BH, but only at the highest 
Kerr order of 4PN. However, we now know that this order is more divergent than lower orders and thus the study 
may not be truly reflective of just how well we can model the Kerr flux function. 

B. Chebyshev flux approximation. 

To model the Chebyshev flux, we defined the interval between the initial and LSO velocities [xi n i,xi so \, where we 
now further define two variables 



X — X[ so T" Xi n i, 
£ Xlso Xi n i- 



From Eqn (|16|) . we have the new iterative relation 



T*(z) = 2 ( T*_ 1 (a;) - T*_ 2 (x), 



(30) 
(31) 



(32) 
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with 



T *(aO = l 



T*{x) 



(33) 



For our study we took Xi n i = 0, however we could just as easily take the velocity when the system is first seen in the 
detector. 

Contrary to the Pade approximation, the Chebyshev re-summation is not effected by null coefficients and can be 
applied to both sums appearing in the factorized flux function. We therefore define the Chebyshev flux function as [l8[ 



F Cn (x) = F N (x) (l-—) 



k=Q 



with the Chebyshev series 



k=0 



k=& 



and 



j2Wk)n(x) = j2f^ 



(34) 



(35) 



(36) 



fc=0 



fe=0 



As the inverted Pade flux was more convergent when compared with the numerical flux, we thus similarly define an 
inverted Chebyshev flux function 



Fwjx) = F N (x) 



•Epole 



fc=0 



(37) 



with 



J2v k (d k )Tt(x) = J2d k a 



k=0 



fc=0 



(38) 



The Chebyshev re-summation was studied in the Schwarzschild case [X 81 ] . and in this case it has already shown a 
more stable behavior which assured the minimization of the error between the approximated and numerical fluxes 
when the order of approximation increased. In general the Chebyshev approximation provided better accuracy at 
lower orders than either the PN or Pade fluxes. The only approximation at which it did not perform the best was 
at the 5.5PN order where it was better than the PN flux but worse than the Pade flux, which for unknown reasons 
is particularly good. We believe that as the Pade approximation had a highly oscillatory nature of performance, 
it just happened to match the numerical flux closely at this PN order. We do not believe that this was a sign of 
convergence of the Pade approximation as the 6PN approximation (if we had it) could again have been worse than 
the PN approximation. While it is not truly reflective as a metric due to the fact that a bad approximation could 
luckily coincide with the numerical flux, the error at the LSO was also investigated for each approximation. While 
the PN and Pade fluxes oscillated wildly in their performance, the Chebyshev approximation essentially reached the 
same constant error value at all orders of approximation which was at least a factor of two better, and sometimes a 
factor of 13 better, than either the PN or Pade flux models. 



IV. RESULTS. 



We know that depending on the value of the spin, the position of the LSO is different. For prograde orbits, the 
test-particle approaches closer to the central black and the system is more relativistic, whereas for retrograde orbits 
the small body falls into the black hole earlier, and the system is less relativistic. Based on this physical consideration, 
we decided to organize our study into three parts : the non-spinning case, the retrograde case and the prograde case. 
In each case, we found the best factorized flux for each re-summation technique, and using these best-choice fluxes 
we compared the two re-summed fluxes with the PN approximation at different orders and for different spins. 
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FIG. 1: Relative error between the approximated fluxes, Fa, and the numerical flux, Fn, in the non-spinning case at 5.5PN for 
the Pade (left panel) and Chebyshev (right panel) flux functions. The dark (black) lines correspond to the factorized flux and 
the light (orange) line correspond to the inverted factorized flux. The solid lines represent a normalization of the logarithmic 
term using xi so and the dashed lines represent a normalization using x pr . The best choice in both cases is the inverted flux 
using xi so for the normalization. 



A. The Schwarzschild Case. 



Our first study concerns the non-spinning case, and before working on any resummation methods we first need to 
find the best factorized flux which will form the basis of these methods. We saw in Section fill Al that two aspects 
of the factorized flux can been tuned : the choice of the normalization in the logarithmic term (we tried xi so and 
x pr ) and the choice to use an inverted or an non- inverted factorized flux. We tried these four cases at 5.5PN and we 
applied on each of them the Pade and the Chebyshev re-summations, the results are shown in Fig[T] The first thing 
that we see is that for both Pade and Chebyshev fluxes, and for the two normalizations, the inverted flux models 
have a smaller error than the non-inverted flux. Depending on which normalization factor is used, the improvement 
ranges between a factor of two and an order of magnitude. This validates the proposition made in Ref 17] to use the 
inverted flux. We also see that for both the Pade and Chebyshev fluxes, using xi so for the normalization, which was 
the choice made in 16], produces a smaller maximum error. Consequently, the best flux, in this sense, is the inverted 
flux using xi so for the normalization. While the results presented here are for the 5.5PN case, we have confirmed a 
similar behaviour at lower orders of approximation. 

Using this factorized flux we compare the PN (black solid curves), Pade (red dot-dashed curves), Chebyshev (green 
doublc-dot-dashed curves) and inverted Chebyshev (blue dashed curves) fluxes with a numerical flux. This comparison 
ranges from 2PN to 5.5PN and is shown in Fig [2] The first thing that we see is that, except at the 2PN order where 
the PN flux has a fortunate coincidence with the numerical flux and the orders where the Pade flux has a singularity, 
the three re-summations are always more convergent than the PN flux. If we compare the performance of each flux 
model as we approach the last stable orbit (and again keeping in mind that this is not on its own the best metric 
of convergence) we can see that the PN approximation returns a relative error e = |1 — Fa/Fn\ (where Fa is the 
approximated flux, and Fn is the numerical flux) of 0.03 < e < 0.4 at various orders of approximation. The Pade 
approximation, when not suffering from a singularity, has relative errors of between 3 x 10~ 2 < e < 8 x 10~ 2 , except 
at the 5.5PN order where it has an error of 1.7 x 10 -3 . The Chebyshev approximation, as expected, tries to minimize 
the maximum error and produce an oscillating error curve. As was shown in Ref [l8j], the Chebyshev approximation 
converges to an almost constant error value of 1.8 x 10 -2 at the LSO. The inverted Chebyshev approximation performs 
fantastically in this case achieving a constant error value of approximately e = 10 -3 at all orders of approximation. 
Therefore, we can see from the figure that the inverted Chebyshev flux at 2PN order is as accurate as the 5.5PN Pade 
approximation (that is difficult to explain) and at least an order of magnitude better than all other models at each 
order of approximation. 

At wider separations, the Chebyshev approximation allows the error to float a little, but in this domain, we do not 
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FIG. 2: Relative error between the approximated flux, Fa, and the numerical flux, Fn, in the non-spinning case, for PN flux 
(black solid curve), the inverted Pade flux (red dash-dotted curve), the Chebyshev flux (green double-dot-dashed curve) and 
the inverted Chebyshev flux (blue dashed curve), for all PN orders from 2PN to 5.5PN. 



require the same level of accuracy as we do closer to the central BH. 



B. Retrograde Orbits. 

As in the Schwarzschild case, we start by searching for the best factorized flux for each re-summation method. We 
demonstrate in the bottom two cells of Fig[3j the comparison between the invcrted/non-invcrtcd fluxes, using different 
normalization constants, at 5.5PN for a spin q — —0.5. The result from this case is representative of the results for all 
other retrograde spins. As in the Schwarzschild case, we see that once again for both the Pade and Chebyshev flux 
models, the inverted flux is the more convergent model. One could argue in the Pade case that the standard flux looks 
to be the best choice as it has a smaller error as we approach the LSO, but it is not the case due to the presence of a 
singularity in the flux. This pole is also present at all other spins and other differing orders of approximation, which 
is why we reject using this flux. We also can see that using x pr is a marginally better choice than using xi so as a 
normalizing velocity. However, the results are so close that one could continue to use the LSO value as a normalizing 
factor. 

Now that we have a model for the optimal factorized flux, we now compare the PN flux (black solid curves), the 
inverted Pade (red dot-dashed curves) and the inverted Chebyshev (blue dashed curves), from 2PN to 5.5PN for the 
four retrograde spin values of q = —0.25, —0.50, —0.75, —0.95, which are shown in Fig HE] The first thing that we 
notice is that the addition of the Schwarzschild and absorption terms to the PN flux to 5.5PN order provides a more 
convergent PN flux model than the standard PN Kerr model at 4PN. This result has the knock-on effect that the 
ensuing 5.5PN factorized flux is better performing than the corresponding factorized flux at 4PN order, thus ensuring 
that the re-summed fluxes starting at 5.5PN order are more convergent than previously studied. 

If we first investigate the performance of the PN flux, we see that this approximation has a slightly improved 
convergence at 5.5PN order as we move from a spin of -0.25 to a spin of -0.95. This is due to the fact that the higher 
retrograde spins being less relativistic and the test-particle plunging further out from the central BH. In general, at 
q = —0.25 the error at the LSO is between 3 x 10 -2 < e < 0.3, improving to between 10~ 2 < e < 0.2 at a spin of 
q = —0.95. While the performance of the PN flux illustrates the oscillatory convergence that we would expect, we 
should point out that as we move towards a maximally counter-rotating spin, the 4.5PN order PN approximation 
performs better and better with an errors of 3 x 10~ 2 , 10~ 2 , 2 x 10~ 4 and 3 x 10 -3 at increasing orders of spin. This 
mimics the unexplainable behaviour we observed in the Schwarzschild case for the Pade approximation at 5.5PN order, 
and in this case the PN flux just happens to have a fortunate coincidence with the numerical flux as we approach the 
LSO. 

The inverted Pade flux does not perform particularly well in the case of retrograde motion. At all spin values, 
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FIG. 3: Relative error between the approximated flux, Fa, and the numerical flux, Fn, in the spinning case at 5.5PN for the 
Pade (left panels) and Chebyshev (right panels) flux functions for spins of q — 0.50 (top panels) and q — —0.50 (bottom panels). 
The dark (black) lines correspond to the factorized flux and the light (orange) lines correspond to the inverted factorized flux. 
The solid lines represent a normalization of the logarithmic term useing xi so , while the dashed lines represent a normalization 
by x pr . For the prograde case the best choices are the the non- inverted flux for the Pade approximation and the inverted flux 
for the Chebyshev approximation, with both fluxes using xi ao for the normalization. For the retrograde case, the best choice is 
the inverted flux using x pr for the normalization, for both Pade and Chebyshev approximations. 



when not suffering from singularities, there are levels of approximation where the Pade approximation is inferior to 
the original PN flux function. This is due to the fact that the factorized flux is divergent at these orders and performs 
much worse than the original PN flux function. In the cases where it does provide an improvement, the overall gain is 
a factor approximately 1.5 over the PN approximation with errors of between 2 x 10~ 2 < e < 0.2 and 10~ 2 < e < 0.1 
at spins of -0.25 and -0.95 respectively. 

The inverted Chebyshev approximation performs consistently well at all spin values, with relative errors of e = 
7 x 10~ 3 at q = —0.25 and e = 10~ 2 at a spin of q = —0.95. Regardless of the spin value, the inverted Chebyshev flux 
is approximately an order of magnitude better at low orders of approximation and at least a factor of two better at 
higher orders. More importantly is the fact that this approximation is almost always better than both the PN and 
inverted Pade fluxes, and essentially converges to the same order of magnitude error at the LSO for all retrograde spin 
values. In the cases where it is less convergent than the other approximations, we again highlight the fact at times 
the PN approximation has a fortunate coincidence with the numerical flux as we approach the LSO, and the Pade 
approximation at 5PN order seems to be very good for high retrograde spin values. However, neither of these facts 
can be considered as global trends, so it is clear that for retrograde orbits, the inverted Chebyshev flux provides a 
more convergent flux model than cither of the other two methods. Furthermore, the inverted Chebyshev flux function 
at 2.5 or 3PN order appears to have a sufficiently small error at all values of separation to be sufficient as a basis for 
a template. 



C. Prograde Orbits. 



We finish our study by looking at the prograde case. As previously, we first search for the best factorized flux. In 
the top two panels of Fig[3]we plot the comparison between the Pade and Chebyshev fluxes at 5.5PN order for q = 0.5. 
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FIG. 4: Relative error between the approximated flux, Fa, and the numerical flux, Fn, with q — —0.25, for the PN (black solid 
curve), inverted Pade (red dash-dotted curve) and inverted Chebyshev (blue dashed curve) flux models at all PN orders from 
2PN to 5.5PN. 
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FIG. 5: Relative error between the approximated flux, Fa, and the numerical flux, Fn, with q — —0.50, for the PN (black solid 
curve), inverted Pade (red dash-dotted curve) and inverted Chebyshev (blue dashed curve) flux models at all PN orders from 
2PN to 5.5PN. 



We see that in both cases the normalization using xi so is marginally better than a normalization using x pr . In the 
Pade case the best choice of flux is the standard Pade flux due to the presence of the poles in both of the inverted flux 
models. On further investigation, these poles appear for all values of prograde spin, so we choose the non-inverted 
factorized flux as the basis for the Pade approximation. On the other hand, for the Chebyshev flux, there is not much 
difference between the four flux models until we approach the LSO. In this case that standard Chebyshev flux, with 
either normalization, begins to diverge. We will therefore use the inverted Chebyshev flux with xi so as a normalizing 
velocity as our preferred model. 

As for the retrograde case we compare the PN (black solid curves), Pade (red dot-dashed curves) and the inverted 
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FIG. 6: Relative error between the approximated flux, Fa, and the numerical flux, Fn, with q — —0.75, for the PN (black solid 
curve), inverted Pade (red dash-dotted curve) and inverted Chebyshev (blue dashed curve) flux models at all PN orders from 
2PN to 5.5PN. 
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FIG. 7: Relative error between the approximated flux, Fa, and the numerical flux, Fn, with q — —0.95, for the PN (black solid 
curve), inverted Pade (red dash-dotted curve) and inverted Chebyshev (blue dashed curve) flux models at all PN orders from 
2PN to 5.5PN. 

Chebyshev fluxes (blue dashed curves), from 2PN to 5.5PN and for q = 0.25, 0.50, 0.75, 0.95, in Fig BHD First of all, 
we can see that, contrary to the retrograde case, when the spin increases the convergence of the flux decreases. This 
is explained by the fact that when the spin increases, the test-particle is now orbiting much closer to the central black 
hole and the system becomes much more relativistic. In this case the PN flux breaks down faster. Before discussing 
the results in full, we should mention here that during our investigation, we noticed that the inverted Chebyshev 
flux was not performing as well as we had expected. When delving deeper into this issue, we discovered that the 
factorized flux function at 5.5PN order behaved much worse than the original PN approximation. The Chebyshev 
approximation, while being able to improve on this somewhat, was unable to recover from the high divergence of the 
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factorized flux. We then studied the behaviour of the factorized flux at different PN orders and discovered that the 
factorized flux at 5PN order was better performing than the 5.5PN PN flux. As the performance of the Chebyshev 
approximation is highly dependent on the starting information, we decided to construct the inverted Chebyshev flux 
using the factorized flux at 5PN order. In the following figures, a new curve represented this truncated inverted 
Chebyshev flux is represented (green double-dot-dashed curves). 

The case of progradc orbits is sufficiently difficult to describe that we divide it into two regimes. We noticed that 
even up to a spin of q — 0.5 the orbits are only slightly more relativistic than the Schwarzschild case, while from 
spins of q = 0.75 onwards, the orbits become truly relativistic to the point that the PN flux begins to break down 
sooner and sooner. However, one global observation is that, even in the prograde case, the PN flux at 4.5 to 5.5PN 
order is more convergent than the standard Kerr 4PN flux. This again justifies the inclusion of the Schwarzschild and 
absorption terms at higher PN orders. 

Firstly, we will discuss the results at spins of 0.25 and 0.5. In both these cases, the PN flux performs as expected 
with fractional errors of between 8 x 10~ 3 < e < 0.5 and 2 x 10~ 2 < e < 0.9 respectively. In the lower spin case the 
error of e = 8 x 10 -3 is due again to a lucky coincidence with the numerical flux as we approach the LSO. In general 
the mean error is on the order of e = 0.1. In the q — 0.5 case, things begin to become a little clearer and we start to 
oscillate wildly between highly convergent and highly divergent. While the Pade flux presents a large improvement on 
the PN flux, in terms of stability, we still see some of the wild fluctuations associated with the PN flux, with relative 
errors of between 4 x 10 -3 < e < 5 x 10~ 2 and 3 x 10 -4 < e < 6 x 10~ 2 . Here as in the PN case, the extremely 
convergent values are due coincidences with the numerical fluxes as we approach the LSO. In general the errors are 
of the order of e = 5 x 10~ 2 and e = 3 x 10~ 2 respectively for both spin values. At this level of spin, the Pade flux 
suffers from singularities at the 3.5PN order which rule out using these orders for template creation. The inverted 
Chebyshev flux starting from the divergent factorized flux at 5.5PN order is almost always better than the PN flux at 
a spin of q = 0.25 achieving an almost constant error of e = 2 x 10~ 2 at the LSO. For this low spin, it in general also 
outperforms the Pade approximation. However for a spin of q = 0.5, while still outperforming the PN flux function, 
the 5.5PN inverted Chebyshev flux is constantly worse than the Pade approximation, achieving a maximum spin error 
of only e = 0.1. This is also reflected at higher spin values and demonstrates the increasing divergence of the 5.5PN 
factorized flux. However, if we now concentrate on the 5PN inverted Chebyshev flux, we can see that it approaches 
an almost constant error value of e = 10~ 2 for a spin of q = 0.25 and e = 5 x 10~ 3 for a spin of q = 0.5. For both 
spin values, it is consistently better than all other approximation models. While in one or two cases it is beaten by a 
fortunate PN or Pade flux model, it clearly provides a more stable model, even at higher values of r/M . 

We next move onto the more relativistic spins of q = 0.75 and q = 0.95. For these two spin values, there is 
not much point in discussing the relative errors for the various approximations as none of them can cope with the 
high relativistic nature of the orbits. In this case, a more indicative metric may be the separation at which the 
relative error reaches e = 0.1. With this in mind, the PN flux is the least performing approximation achieving the 
threshold error at separations of between 5 < r/M < 10 for both spins. In each case the PN flux at 5.5PN performs 
exceptionally well, getting as close as r — 3.5M before having an error of 10 _1 . However, this is clearly not reflective 
of convergence, as a 6PN flux could perform much worse. The Pade flux again improves over the PN flux, except at 
the 5.5PN order, reaching the threshold error at separations of between 3<r/M<5in both cases, with the odd 
exception performance in the highest spin case at 5PN where a separation of r = 2.5AT is reached. However, the 
Pade flux approximation now begins to suffer from singularities not only at 3.5PN order as we saw for lower spins, 
but also now at the 3PN order. The 5.5PN inverted Chebyshev flux has a relative error of 10 _1 at separations of 
around r = 5M in both cases. While it achieves a consistent separation before reaching the threshold, it can only be 
marginally considered to be better than the PN flux, and is constantly outperformed by the Pade flux. However, the 
5PN inverted Chebyshev flux reaches the threshold value at separations of approximately r = AM in both spin cases. 
While sometimes outperformed at a particular approximation by either the PN or Pade flux models, we believe that 
the 5PN inverted Chebyshev flux represents the most stable and consistently performing model for prograde orbits. 



V. CONCLUSION. 

In this work we extended the studies of [H-18] by taking into account the absorption of the GW by the central 
black hole, and by adding the information from the non-spinning case, we created a 5.5PN Kerr flux. We saw that, 
because the 4PN flux is divergent, adding these higher order terms improves the convergence of the PN flux in all 
cases. 

By first looking at the non-spinning case, we searched for the most convergent factorized flux at 5.5PN for each 
re-summation method. We tried a number of different combinations : the standard factorized flux given in Ref. fl6j 
or the inverted factorized flux proposed in Ref. [l7| . and we also tried xi so , the invariant speed at the last stable 
orbit and x pr , the invariant speed at the photon ring, as a normalization factor in the logarithmic term of the flux. 
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FIG. 8: Relative error between the approximated flux, Fa, and the numerical flux, Fm, with q = 0.25, for PN (black solid 
curve), Pade (red dash-dotted curve), inverted 5.5PN Chebyshev (blue dashed curve) and inverted 5PN Chebyshev (green 
double-dot-dashed curve) flux models at all PN orders from 2PN to 5.5PN. In the last panel, the green double-dot-dashed curve 
is again the inverted 5PN Chebyshev flux as a comparison. 
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FIG. 9: Relative error between the approximated flux, Fa, and the numerical flux, Fm, with q = 0.50, for PN (black solid 
curve), Pade (red dash-dotted curve), inverted 5.5PN Chebyshev (blue dashed curve) and inverted 5PN Chebyshev (green 
double-dot-dashed curve) flux models at all PN orders from 2PN to 5.5PN. In the last panel, the green double-dot-dashed curve 
is again the inverted 5PN Chebyshev flux as a comparison. 

It appeared that for both Pade and Chebyshev flux models, the best choice was to use the inverted flux normalized 
using xi so , which confirms the choices made in [l6| and [13], respectively. Then we compared the two resummation 
methods and the PN flux, at different orders, with a numerical flux. We observed that, except when it presents a 
pole, the Pade approximation is in general better than the PN flux. However, while the standard Chebyshev flux 
performed extremely well, the better approximation was clearly the inverted Chebyshev flux function. 

Next we studied the case of retrograde orbits, and once again we started by searching for the best factorized flux. 



15 




FIG. 10: Relative error between the approximated flux, Fa, and the numerical flux, Fm, with q — 0.75, for PN (black solid 
curve), Pade (red dash-dotted curve), inverted 5.5PN Chebyshev (blue dashed curve) and inverted 5PN Chebyshev (green 
double-dot-dashed curve) flux models at all PN orders from 2PN to 5.5PN. In the last panel, the green double-dot-dashed curve 
is again the inverted 5PN Chebyshev flux as a comparison. 
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FIG. 11: Relative error between the approximated flux, Fa, and the numerical flux, Fm, with q — 0.95, for PN (black solid 
curve), Pade (red dash-dotted curve), inverted 5.5PN Chebyshev (blue dashed curve) and inverted 5PN Chebyshev (green 
double-dot-dashed curve) flux models at all PN orders from 2PN to 5.5PN. In the last panel, the green double-dot-dashed curve 
is again the inverted 5PN Chebyshev flux as a comparison. 



Contrary to the previous case, for negative spins the marginally better normalization choice is using x pr , and as in 
the Schwarzschild case, the inverted flux is the best choice for the two techniques. When comparing the different 
fluxes we saw that, except for a few particular cases, the re-summation methods are better than the PN flux. At 
times however, the PN flux has a fortunate crossing with the numerical flux near the LSO and therefore appears 
to work well. However, the fluctuations is convergence from one order to another means that it is difficult to rely 
on the PN approximation. In general, the Pade approximation improves on the PN flux, but again there are orders 
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of approximation where it performs worse than the PN approximation. The inverted Chebyshev approximation 
outperforms both other flux models, improving the accuracy of the flux model by a factor of two. 

Our final study involved prograde orbits and once again we found that the best normalization is using Xi so , but 
while the inverted flux was the best choice for the Chebyshev approximation, it appeared that the standard flux is a 
better choice for the Pade flux. Here we found that the PN approximation performs worse and worse as we increase the 
spin magnitude, except at the 5.5PN order where it performs surprisingly well. The Pade flux is a great improvement 
over the PN flux, but again oscillates between good and bad performance. In this case, we found that the inverted 
Chebyshev flux at 5.5PN order did not perform as well as expected. Upon further investigation, we found that this 
was due to the factorized flux being divergent at this order. However, by using the 5PN factorized flux, we were once 
again able to construct a very convergent Chebyshev flux model. While more marginal, the inverted Chebyshev flux 
is the most stable model, converging to the same order of error at the LSO. 

This work improves upon previous studies and introduces a new inverted Chebyshev flux model. This model is 
highly convergent and stable. In most cases, this flux model is more convergent at the 2.5PN level than either the 
PN or Pade fluxes at 5.5PN order. We believe that this model provides a good basis for the creation of theoretical 
templates for GW astronomy. 
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Appendix A: PN flux coefficients. 



Here we present the post-Newtonian coefficients taking into account the non-spinning terms up to x 11 , the spinning 

terms up to a; 8 and the absorption terms up to x . 
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where 7 is the Euler's number, k = \J\ — q 2 , and the A n , B n , and C n are defined by 
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where (z) and ip^ (z) are respectively the digamma and the trigamma functions. For ease and speed of calculation, 
the digamma and trigamma functions were generated using a Chebyshev-Pade approximation such that the maximum 
error between the approximation and the analytic expressions is less than 10~ 10 for all spin values. 



